# 0.0 Load libraries ----
packages<-c('tidyverse','multcomp','rio')
lapply(packages,require, character.only=TRUE)
# 1.0 Load and format data ----
data_path = ''
ds0 = import(data_path)
ds = ds0 %>% 
  mutate(ew_treatment_m = factor(ew_treatment_m, levels=c("control", "creditclaim", "porkbusting")),
         issuematch_f = relevel(factor(issuematch_f), ref="control"))
# 2.0 Hypothesis 1 Test ----
# 2.1 H1 Visual Mean Comparison ----
ds %>% 
  drop_na(ew_therm, ew_treatment_m) %>% 
  group_by(ew_treatment_m) %>%
  summarise(nrow = length(ew_therm), 
            mean = mean(ew_therm), 
            sd = sd(ew_therm), 
            se = sd/sqrt(nrow), 
            ci = 1.96*se) %>% 
  ggplot(., aes(x = ew_treatment_m, y = mean, shape = ew_treatment_m)) + 
  geom_pointrange(aes(ymin = mean - ci, ymax = mean+ci), size = 1) + 
  scale_shape_manual(values=c(19, 15, 17), labels = c("Control", "Credit Claiming", "Pork Busting"))+
  scale_y_continuous(limits = c(55,70)) + 
  labs(title = "(A) Avg. Approval (0-100)",
       y = "",
       x = "")
# 2.2 H1 ANOVA w/ Tukey HSD ----
h1anova = aov(ew_therm ~ ew_treatment_m, ds)
mod3.1 = glht(h1anova, linfct=mcp(ew_treatment_m = "Tukey"))
mod3.2 = summary(mod3.1)$test
mod3.2lwr = format(round(as.data.frame(confint(mod3.1)$confint)$lwr, 3), nsmall=3)
mod3.2upr = format(round(as.data.frame(confint(mod3.1)$confint)$upr, 3), nsmall=3)
mod3.2ci = paste0("[", mod3.2lwr, ", ", mod3.2upr,"]")
mod3.3 <- cbind(round(mod3.2$coefficients, 3), round(mod3.2$sigma, 3), mod3.2ci, format(round(mod3.2$pvalues, 3), nsmall = 3))
colnames(mod3.3) <- c("Estimate", "Std. Error", "95% Confidence Interval", "p-value")
mod3.3
# 3.0 Hypothesis 2 Test ----
# 3.1 H2 Visual Mean Comparison ----
ds %>% 
  drop_na(ew_eff, ew_treatment_m) %>% 
  group_by(ew_treatment_m) %>%
  summarise(nrow = length(ew_eff), 
            mean = mean(ew_eff), 
            sd = sd(ew_eff), 
            se = sd/sqrt(nrow), 
            ci = 1.96*se) %>% 
  ggplot(., aes(x = ew_treatment_m, y = mean, shape=ew_treatment_m)) + 
  geom_pointrange(aes(ymin = mean - ci, ymax = mean+ci), size=1) + 
  scale_shape_manual(values=c(19, 15, 17), labels = c("Control", "Credit Claiming", "Pork Busting"))+
  scale_y_continuous(limits = c(5.5,7.0)) + 
  labs(title = "(B) Avg. Lawmaking Effectiveness (0-10)",
       y = "",
       x = "") 
# 3.2 H2 ANOVA w/ Tukey HSD ----
h2anova = aov(ew_eff ~ ew_treatment_m, ds)
mod4.2 = glht(h2anova, linfct=mcp(ew_treatment_m = "Tukey"))
mod4.2.2 = summary(mod4.2)$test
mod4.2lwr = format(round(as.data.frame(confint(mod4.2)$confint)$lwr, 3), nsmall=3)
mod4.2upr = format(round(as.data.frame(confint(mod4.2)$confint)$upr, 3), nsmall = 3)
mod4.2ci = paste0("[", mod4.2lwr, ", ", mod4.2upr,"]")
mod4.3 <- cbind(format(round(mod4.2.2$coefficients, 3), nsmall=3), round(mod4.2.2$sigma, 3), mod4.2ci, format(round(mod4.2.2$pvalues, 3), nsmall = 3))
colnames(mod4.3) <- c("Estimate", "Std. Error", "95% Confidence Interval", "p-value")
mod4.3
# 4.0 Hypothesis 3 Test ----
# 4.1 H3 Visual Mean Comparison ----
ds %>% 
  drop_na(ew_resp, ew_treatment_m) %>% 
  group_by(ew_treatment_m) %>%
  summarise(nrow = length(ew_resp), 
            mean = mean(ew_resp), 
            sd = sd(ew_resp), 
            se = sd/sqrt(nrow), 
            ci = se*1.96) %>% 
  ggplot(., aes(x = ew_treatment_m, y = mean, shape=ew_treatment_m)) + 
  geom_pointrange(aes(ymin = mean - ci, ymax = mean+ci), size=1) + 
  scale_shape_manual(values=c(19, 15, 17), labels = c("Control", "Credit Claiming", "Pork Busting"))+
  scale_y_continuous(limits = c(5.5,7.0)) + 
  labs(title = "(C) Avg. Fiscal Responsibility (0-10)",
       y = "",
       x = "")
# 4.2 H3 ANOVA w/ Tukey HSD ----
h3anova = aov(ew_resp ~ ew_treatment_m, ds)
mod5.2 = glht(h3anova, linfct=mcp(ew_treatment_m = "Tukey"))
mod5.2.2 = summary(mod5.2)$test
mod5.2lwr = format(round(as.data.frame(confint(mod5.2)$confint)$lwr, 3), nsmall=3)
mod5.2upr = format(round(as.data.frame(confint(mod5.2)$confint)$upr, 3), nsmall=3)
mod5.2ci = paste0("[", mod5.2lwr, ", ", mod5.2upr,"]")
mod5.2.3 <- cbind(format(round(mod5.2.2$coefficients, 3), nsmall=3), round(mod5.2.2$sigma, 3), mod5.2ci, format(round(mod5.2.2$pvalues, 3), nsmall = 3))
colnames(mod5.2.3) <- c("Estimate", "Std. Error", "95% Confidence Interval", "p-value")
mod5.2.3
# 5.0 Hypothesis 4 Test ----
# 5.1 H4 Difference-in-means test ----
ds_c_r = ds %>% filter(partyr=="rep" & ew_treatment_m=="control")
ds_c_d = ds %>% filter(partyr=="dem" & ew_treatment_m=="control")
ds_cc_r = ds %>% filter(partyr=="rep" & ew_treatment_m=="creditclaim")
ds_cc_d = ds %>% filter(partyr=="dem" & ew_treatment_m=="creditclaim")
t.test(ds_cc_r$ew_therm, ds_c_r$ew_therm, alternative = "two.sided", var.equal = FALSE)
t.test(ds_cc_d$ew_therm, ds_c_d$ew_therm, alternative = "two.sided", var.equal = FALSE)
# 5.2 H4 Visualization ----
ds %>% 
  drop_na(partyr, ew_therm) %>%
  filter(ew_treatment_m!="porkbusting") %>% 
  group_by(partyr, ew_treatment_m) %>% 
  summarise(nrow = length(ew_therm), 
            mean = mean(ew_therm), 
            se = sd(ew_therm)/sqrt(nrow),
            ci = 1.96*se) %>% 
  ggplot(., aes(x = partyr, y = mean, shape = ew_treatment_m)) + 
  geom_pointrange(aes(ymin = mean-ci, ymax = mean+ci), position = position_dodge(width = .5), size = 1) + 
  labs(title = "(A) Control-Credit Claiming", 
       shape = "",
       x = "", 
       y = "Avg. Approval (0-100)") + 
  scale_x_discrete(labels = c("Democrats\n(p>0.05)", "Republicans*\n(p<0.05)")) + 
  scale_shape_manual(values=c(19, 15), labels = c("Control", "Credit Claiming"))+
  scale_y_continuous(limits=c(55, 71)) + 
  scale_color_manual(values = c("darkgrey", "black"), guide='none')
# 6.0 Hypothesis 5 Test ----
# 6.1 H5 Difference-in-means test ----
ds_pb_r = ds %>% filter(partyr=="rep" & ew_treatment_m=="porkbusting")
ds_pb_d = ds %>% filter(partyr=="dem" & ew_treatment_m=="porkbusting")
t.test(ds_pb_r$ew_therm, ds_c_r$ew_therm, var.equal = FALSE)
t.test(ds_pb_d$ew_therm, ds_c_d$ew_therm, var.equal = FALSE)
# 6.2 H5 Visualization ----
ds %>% 
  drop_na(partyr, ew_therm) %>%
  filter(ew_treatment_m!="creditclaim") %>% 
  group_by(partyr, ew_treatment_m) %>% 
  summarise(nrow = length(ew_therm), 
            mean = mean(ew_therm),
            se = sd(ew_therm)/sqrt(nrow), 
            ci = se*1.96) %>% 
  ggplot(., aes(x = partyr, y = mean, shape = ew_treatment_m)) + 
  geom_pointrange(aes(ymin = mean - ci, ymax = mean+ci), position = position_dodge(width = .5), size = 1) + 
  labs(title = "(B) Control-Pork Busting", 
       shape = "",
       x = "", 
       y = "") + 
  scale_x_discrete(labels = c("Democrat*\n(p<0.05)", "Republican\n(p>0.05)")) + 
  theme(legend.position = "bottom", legend.background = element_rect(color = NA)) + 
  theme(legend.position="right") + 
  scale_shape_manual(values=c(19, 17), labels = c("Control", "Pork Busting")) +
  scale_y_continuous(limits=c(55, 71)) + 
  scale_color_manual(values = c("black", "darkgrey"), guide='none')
# 7.0 H4/H5 In-party/Out-party Analysis ----
# 7.1 H4/H5 IO difference-in-means test ----
ds_io_c = ds %>% filter(ioMessage=="control")
ds_io_o = ds %>% filter(ioMessage=="outM")
ds_io_i = ds %>% filter(ioMessage=="inM")
t.test(ds_io_c$ew_therm, ds_io_i$ew_therm, var.equal = FALSE)
t.test(ds_io_c$ew_therm, ds_io_o$ew_therm, var.equal = FALSE)
# 7.2 H4/H5 IO Visualization ----
ds %>% 
  drop_na(ioMessage, ew_therm) %>%
  group_by(ioMessage) %>% 
  summarise(nrow = length(ew_therm), 
            mean = mean(ew_therm),
            se = sd(ew_therm)/sqrt(nrow), 
            ci = 1.96*se) %>% 
  ggplot(., aes(x = ioMessage, y = mean)) + 
  geom_pointrange(aes(ymin = mean - ci, ymax = mean+ci), position = position_dodge(width = .5), size = 1, shape =1) + 
  labs(title = "(C) In-party/Out-party Status", 
       shape = "",
       x = "", 
       y = "") + 
  scale_x_discrete(labels = c("Control", "In-party\n(p>0.05)", "Out-party*\n(p<0.05)")) + 
  scale_y_continuous(limits=c(55, 71))
# 8.0 Regression Models ----
# 8.1 Estimate Regression Models ----
summary(lm(ew_therm ~ 
               issuematch_f +
               partyr +
               partymatch +
               gender +
               poly(age,4) +
               white +
               black +
               hisplat +
               asam +
               natam +
               rother + 
               education +
               income, 
             ds%>%drop_na(age)))
summary(lm(ew_eff ~ 
               issuematch_f +
               partyr +
               partymatch +
               gender +
               poly(age,4) +
               white +
               black +
               hisplat +
               asam +
               natam +
               rother + 
               education +
               income, 
             ds%>%drop_na(age)))
summary(lm(ew_resp ~ 
               issuematch_f +
               partyr +
               partymatch +
               gender +
               poly(age,4) +
               white +
               black +
               hisplat +
               asam +
               natam +
               rother + 
               education +
               income,
             ds%>%drop_na(age)))
